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Abstract 

The phenomenon called the "sub-Sun" is the specular reflection of sunlight by hori- 
zontally oriented plates of ice. Although well-known in meteorological optics, the hydrody- 
namics of the orientation is not quantitatively understood. I review the theory of torques 
on objects at low Reynolds numbers, define coefficients C Q , C p , and which describe 
the orienting torques on discs, rods, and hexagonal prisms, and report here the results of 
experiments to measure C a and C p . 
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I. Introduction 

From an airplane it is sometimes possible to see a nearly specular reflection of the Sun 
from a flat layer of stratus clouds, known in meteorological optics (1-3) as a "sub-Sun". 
This striking phenomenon is not unusual; I have seen it twice in a few hundred flights (on 
most of which the lighting or meteorological conditions were unsuitable) in the last six 
years. On the occasions on which I observed a sub-Sun (early morning in March over New 
Mexico, and late afternoon in June over Ohio) it was broadened by less than a degree of 
arc. Where the stratus was interrupted by cumulus the sub-Sun disappeared, and where 
the stratus was partly transparent and a body of water appeared below it, the nearly 
specular reflection from the stratus was superimposed on the sharper specular reflection 
from the water. 

The well-known (1-3) explanation of the sub-Sun is that the stratus clouds consist 
of small thin plates of ice maintaining a horizontal orientation as they fall. The angular 
broadening reflects a small dispersion in their orientations; diffraction may also contribute 
but is typically smaller. If diffraction were dominant the sub-Sun would be larger in red 
than in blue light, appearing blue at the center but surrounded by a red halo, while angular 
dispersion of the plates leads to no dependence on color. I noticed no color dependence. 

There are related phenomena in meteorological optics. Light pillars (4), like the sub- 
Sun, require horizontally oriented ice plates. The better-studied halos (5,6) may be caused 
by oriented or unoriented ice crystals, either plates or columns. The principles of the 
hydrodynamics of the orientation of small falling objects are understood (7,8), but the 
formidable formalism has been neither explicitly evaluated (except for spheroids of small 
eccentricity and long rods) nor tested experimentally. 

Under what conditions will a small, thin falling plate of ice maintain an accurately 
horizontal orientation? At high Reynolds number Re>la falling plate leaves a turbulent 
wake (9,10). Its center of drag lies close to its leading edge or surface; any steady orientation 
is unstable; it tumbles, and its path is irregular because of large horizontal forces arising 
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during its tumbling (ice crystals, unlike airplanes, rockets, and arrows, are not equipped 
with stabilizing tails!). This is readily verified by dropping a penny into a jar of water, 
an experiment in which Re ps 3000; it tumbles and usually hits the sides. For Re ps 100 a 
falling disc may oscillate periodically about a horizontal orientation as it leaves behind a 
regular vortex street. This may be seen by dropping aluminum foil discs of various radii into 
water. At Re < 100, however, tumbling and oscillations are strongly damped by viscosity. 
Intuitive concepts from our everyday experience with high-Re flows are still qualitatively 
applicable, and show that a vertical orientation (edge-on to the flow) is unstable; if the 
plate tilts the hydrodynamic force on its leading edge acts to amplify the tilt. However, 
the horizontal orientation (face-on to the flow) is stable; if the plate tilts the wake of the 
leading edge partly shields the trailing edge from the flow, reducing the drag on it; the 
resulting torque restores the horizontal orientation and the disturbance is quickly damped. 

II. Low Reynolds Number Flow 

Very small plates will fall slowly, with Re 1. The theory of low- Re flow is presented 
by Happel and Brenner (11), who include most of the results used in this section for flow 
in the limit Re — > 0. In the Navier-Stokes equation 

g + (<;.VK=--Vp + ^V 2 <7, (1) 
at p p 

where v is the fluid velocity field, p its density, p the pressure, and rj its dynamic viscosity, 
the nonlinear inertial term (y ■ V)v ~ v 2 /£, where I is a characteristic length, is O(Re) 
times smaller than the terms on the right hand side, and may be neglected. Similarly, 
the time derivative term ^ ~ v 2 /£ is of the same order of smallness (in the absence of 
time-dependent external forcing), and also may be neglected. This leaves the creeping flow 
equation 

V p = v V 2 v, (2) 

which is linear in velocity. Because of the linearity of (2), linear combinations of solutions 
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which satisfy homogeneous boundary conditions (such as v = on a solid boundary) are 
also solutions. 

An ellipsoidal body sinking under gravity in the limit Re — > suffers no hydrodynamic 
torque, regardless of its orientation: Decompose the fluid velocity at infinity into compo- 
nents along the ellipsoid's principal axes. Each of these corresponds to a solution (If) of 
(2) which satisfies a zero- velocity boundary condition on the body's surface, has a uniform 
velocity at infinity, and (by symmetry) exerts no hydrodynamic torque about its geometric 
center. Thus the body suffers no hydrodynamic torque, whatever its orientation, and, by 
the symmetry of the solutions to (2), the center of drag is also its geometric center, ff this 
is also its barycenter (as in the case of an ellipsoid of uniform density), then there is no 
net torque, and all orientations are neutrally stable. 

The drag tensor D of a spheroid or triaxial ellipsoid, defined by the relation between 

— * 

the applied force F and its velocity u 

F = D • u, (3) 

is not isotropic, ff an external force is not directed along a principal axis then u will 
not be parallel to F; in a gravitational field an obliquely oriented ellipsoid will not sink 
straight down. An elementary calculation, using the classic results (11-13) for a circular 
disc F|| = -yryr-uy and F± = 16rjru_\_, where the subscripts indicate velocities and forces 
parallel and perpendicular to the disc's surface, yields the result for the angle cf> between 
the velocity vector and the nadir: 

i / sin 9 cos 9 \ . A . 

= tan" 1 =- , (4) 

V2 + sin 2 67 V ' 

where 9 is the angle between the surface normal and the vertical. The maximum value of 
is tan-^f^ 1 / 2 ) « 11.5°, which is found for 9 = sin" 1 ((2/5) 1 / 2 ) « 39.2°. 

ft is possible to generalize (in the limit Re — > 0) the result of zero torque from ellipsoids 
to a larger class of shapes. The general relation (11,14) between torque r and body velocity 
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u for a nonrotating body is 

r = T-u, (5) 

where T is a second rank pseudotensor (it changes sign under inversion of all coordinates) 
because f is a pseudovector. It is not possible to construct a nonzero second rank pseu- 
dotensor from a nonchiral geometric shape. Hence T = and there is no hydrodynamic 
torque on any nonchiral object in the limit Re — > 0. It is, however, possible for the centers 
of drag and mass to differ, so that the combination of hydrodynamic and gravitational 
forces may lead to a net torque and a preferred orientation; a nail is an example of such 
an object, which is nonetheless nonchiral. 

III. Not Quite Such Low Reynolds Number Flow 

The observations of the sub-Sun clearly require that flat plates of ice maintain a 
horizontal orientation, but we have seen that in the limit Re — > there can be no aligning 
torque. The explanation is that for any finite Re there is an aligning torque, whose 
magnitude may be estimated as first order in Re and in 9: 

re = C F ± r9Re r , (6) 

where C a is a dimensionless coefficient of order unity applicable to discs and flattened 
oblate spheroids. Applying F± for a horizontally falling disc of radius r and half-thickness 
h to a thin hexagonal plate, and using the resulting u = irgrhp/(8r]) and Re r = urjv = 
irgr 2 hp/(8r]v), where v is the kinematic viscosity, yields 

„2fi 2 5/,2 2 

4 r]u 

A torque of magnitude r acting on a disc produces a rotation rate (11,13,15) (about 
any axis!) in low Re flow 9 = 3re / (32i]r 3 ) . Using this result as an approximation for a thin 
hexagonal plate, and combining the rotational and translational flows (permitted in the 
Re — > limit because both satisfy the homogeneous boundary condition of zero relative 
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velocity on the solid body) yields an exponential decay time t of any deviation from the 
horizontal orientation: 

t g 128 

6 3n 2 C g 2 r 2 h 2 p 2 ' {) 

In this time the plate falls a distance d = ut, conveniently expressed in terms of its radius: 

d - 2 

r ~ 3C Re r ' ^ ^ 

A steady shear of the wind velocity produces a tilt of a falling plate 6 w where 
is an appropriate component of the strain rate tensor of the wind shear. In a stratus cloud 
the strain rate is likely to be small or negligible if the air is stably stratified. However, 
convective turbulence can readily disrupt the orientation of the plates, accounting for the 
observation that the specular sub-Sun disappeared wherever the stratus was interrupted 
by cumulus clouds. 

The dominant influence disrupting the alignment of small plates in still air is orienta- 
tional Brownian motion (2,13-16). The angular diffusion coefficient D Q = 3/c_BT/(32?yr 3 ). 
The root-mean-square dispersion of the angular orientation about a single horizontal axis 
is 

For r = 30/U, h = l.bfi and C a = 1 the dispersion (6* 2 ) 1 / 2 0.01 radian, roughly the largest 
value permitted by the observations. For these dimensions t 0.2/C o sec, u 0.9 cm/sec, 
d 0.1 cm ~ 30r, and Re r ss 0.02. Significantly smaller particles would not produce a 
clear specular reflection because of their Brownian angular dispersion, while much larger 
ones would precipitate rapidly. 

Ice crystals may be needle-like, and their behavior is also of interest. The drag forces on 
prolate ellipsoids (11,13,15) approximate those on needle-like hexagonal cylinders. At low 
Reynolds numbers they are Fy = 4irr)ru\\/(\n(2r/h) — \) and F± = 8n:riru±/(ln(2r /h) + 1), 
where r is the longest semi-axis (half the needle's length) and h is each of its shorter semi- 
axes (the needle's radius), only the leading terms in h/r are taken, and the subscripts 
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indicate velocities and forces parallel and perpendicular to the needle's length. An ele- 
mentary calculation of the direction of gravitational settling yields, in analogy to equation 
(4), 



( ^ = tan- 1 ( \ : _ 2 l\ — — ), (11) 



_ x f (1 — C) sin 6* cos 9 
sin 2 6* + ( cos 2 6* 

where 6* is the angle between the longest axis and the horizontal and ( = (ln(2r//i) + 
|)/2(ln(2r//i) — |). For C = § equation (4) is recovered; plausible values of ( are quite 
close to | and lead to a very similar dependence of 4> on 

In the limit Re — > needles will fall with indeterminate orientation, by the same 
argument which applies to flat plates. At small but finite Re there will similarly be an 
aligning torque tending to make their long axes horizontal 

Te = C p F ± r6Re r , (12) 

where the coefficient C p applies to needles and slender prolate spheroids. The results 
analogous to equations (7)-(10) are 

Te = 2 i^ - ' (13) 



16 rj 2 v 



3C p (\n(2r/h) - \){\n{2r/h) + \) g 2 h 4 p 2 ' 

d _ ]n(2r/h) + \ 

r ~ 3C p (ln(2r/h) - |)Re r ' 



(14) 
(15) 



1^1/2 = ( 4 k B T W \ 1/2 

{ ; UC p (ln(2r//i) + i) (/VMp 2 ,/ ' 1 J 

and the angular diffusion coefficient of prolate spheroids D p = 3(ln(2r//i) — ^)ksT/ (87r?yr 3 ). 

A needle-like ice crystal may also rotate about its long axis. Because it is a hexagonal 
prism, there are preferred values of its orientation angle ip about this axis. In analogy to 
the argument for thin plates, the preferred orientation is that in which two prism faces 
(top and bottom) are horizontal if 9 = 0. The previous discussion may also be applied 
to this rotation. Because a hexagonal prism is close to a circular cylinder, for its rotation 
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about its long axis I use the relation (12) for a long cylinder of radius h and length 2r: 
ip = T^j '{9>Txrh 2 rf) (rather than that for an ellipsoid). The torque is given by 

= C^F±htfjRe h , (17) 

where h is used in place of r both in the lever arm and in the Reynolds number Re^; the 
coefficient is expected to be small because of the resemblance of a hexagonal prism to 
a circular cylinder, for which = 0. The results analogous to equations (7)-(10) and 
(13)-(16) are 

7rCV(ln(2r//Q + l) , ia * 
T * = 2 (18) 



16 rfu 



C^(ln(2r//i) + i)^/iV' 
d _ ln(2r//i) + \ 



h C^Reh 



(19) 
(20) 



W; U^(ln(2r//i) + i)<7 2 r/,V; ' 1 J 

and the angular diffusion coefficient = ksT / {Sitrirh 2 ). 

For an ice needle with r = 30/i, h = 3/j, and C p = 1, (6* 2 ) 1 / 2 ~ 0.04 radian, while 

(ip 2 ) 1 / 2 ps 0.4/Cy 2 radian. Brownian angular dispersion is larger for needles than for 

plates of comparable dimensions, and is particularly disruptive of the orientation of needles 

about their long axes. However, a slender needle is lighter and falls more slowly than a 

plate whose diameter equals the needle's length, so longer needles may be found in clouds, 

and these longer needles are better oriented. 

IV. Experiments 

The coefficients C G , C p and can be measured, either in the laboratory or from 
numerical experiments with three-dimensional hydrodynamics calculations. I therefore 
undertook experimental measurements. These turned out to be unexpectedly difficult, in 
part because the coefficients are small. For Re in the range 1-100 the orienting torques 



on falling discs and rods are large, and they become horizontal very quickly (in a few 
diameters). However, C a , C p and are defined as coefficients in a low- Re expansion of 
the torque, so that experiments must be done for Re 1, for which the orienting torques 
are small. This makes the experiments sensitive to small imperfections and heterogeneities 
in the shape and mass distributions of the test objects. The requirement that the orienting 
torques be large enough to observe (and be large compared to spurious torques resulting 
from imperfections) conflicts with the requirement that the first term in an expansion in 
Re be a good approximation. In the experiments reported here 0.10 < Re r < 0.25, with 
most data obtained at the upper end of that range. 

It is easy to obtain very low Re by using very viscous fluids, but the orienting torques 
would be too small to observe. For larger Re high viscosities imply large test objects (and 
increased wall effects) or higher speeds (making observations more difficult). Obtaining 
the required Re with test bodies of reasonable size requires a quite modest viscosity; I used 
soybean oil, for which r] = 0.6 poise at room temperature (17). The size of test bodies also 
represents a compromise; very small objects are difficult to make, but large ones would 
require an excessively large fluid volume to avoid wall effects on the flow, which decrease 
linearly in the ratio of body size to wall distance (11); I used discs with radii in the range 
0.16-0.32 cm in a full tank of dimensions 15 cm x 20 cm x 30 cm. The modest viscosity 
requires very light test bodies to obtain Re in the desired range; the discs were punched 
using ordinary paper punches from hard 1 mil (0.0025 cm) thick aluminum foil (household 
aluminum foil is softer, and discs punched from it have more ragged edges). Data were 
recorded using a hand-held 35 mm camera at about 50 cm distance. 

When working with small and light bodies, surface forces are important. Attempts 
to use sugar-water solutions as the fluid failed because air bubbles invariably attached 
themselves to the test bodies when they were pushed through the air-fluid surface. This 
problem was not solved by the addition of surfactants, and required the use of oil, which 
better wets aluminum and pencil leads. Even so, small bubbles of air often had to be 
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dislodged by stirring the fluid near them after the test bodies were immersed. 

Data were obtained from sinking discs with r = 0.32 cm (j" diameter), for which 
Re r = 0.23, and from discs with r = 0.24 cm diameter), for which Re r = 0.13. 

Measurements of eight discs with r = 0.32 cm and of six discs with r = 0.24 cm led to the 
results 

/ 0.169 ±0.017 Re r = 0.23 . , 

°° ~ \ 0.26 ± 0.05 Re r = 0.13. [ > 

The (geometric) mean 9 during the experiments was 16° for the r = 0.32 cm discs and 

23° for the r = 0.24 cm discs; Eq. (22) includes a finite 9 correction factor 26>(cos 2 9 + 

| sin 2 9) I sin 29 which corrects for the dependence of the vertical component of u on 9 and 

also replaces 9 in Eq. 6 by \ sin 29, giving the correct limiting torques at 9 = 90° (r = 0) 

and at 9 = 0. 

The experiments with r = 0.32 cm also yielded three discs whose relaxation to hori- 
zontal was strikingly slower (by about 50%) , than that of the eight included above (whose 
dispersion was about 10%). These outliers could have been the consequence of unobserved 
small attached bubbles or manufacturing imperfections (punching leaves a region of plastic 
flow and tearing around a disc's edge which is not exactly symmetric, and is sometimes 
grossly asymmetric) which displace the center of gravity from the center of drag. If these 
outlying data points are included the result is C a = 0.154 ±0.031, not significantly different 
than the result above. In the experiments with r = 0.24 cm two discs did not relax to the 
horizontal at all and were rejected entirely; the dispersion of relaxation rates among the 
remainder was greater than for the larger discs, but there were no remaining outliers. 

The errors quoted are the standard deviations of the measurements of individual 
discs. If the scatter were Gaussian with zero mean, then the uncertainty in C Q would 
be about 0.4 of that quoted, but this assumption is probably unduly optimistic. The 
measurements at two different Reynolds numbers appear to disagree significantly; known 
sources of systematic error are believed to be much less than the quoted errors. This may 
imply that the Re — > limit is not adequately approximated by Reynolds numbers of 0.2. 
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It would be desirable to extend the measurements to yet lower Reynolds numbers, but 
smaller discs are both more sensitive to imperfections (because of their lower Reynolds 
numbers and smaller aligning torques) and harder to observe accurately, as is shown by 
the greater dispersion of the results for r = 0.24 cm discs. 

Cylinders of mechanical pencil "lead" with h = 0.0183 cm (sold as 0.3 mm diameter) 
and lengths 0.44-0.66 cm were used to determine C p . Despite the apparent uniformity of 
the cylindrical shape (their ends were flattened by mounting them in a jeweler's lathe and 
grinding them against a flat grindstone) the results were rather scattered. The origin of 
the scatter may be heterogeneity in density or remaining imperfections of the end faces. 
The result was 

C p = 0.059 ± 0.016, (23) 

at a mean Re r = 0.16. As above, the quoted error is the dispersion of the (seven) individual 
measurements only. This rough result is less than the theoretical prediction of C p ~ 
0.12 for the rods used (18,19), and again may reflect either manufacturing imperfections, 
unobserved attached bubbles, or a failure to reach the Re — > regime at finite Re. 
is much more difficult to measure, requiring precision hexagonal cylinders, and I did not 
attempt it. 

I thank D. R. Nelson and H. A. Stone for discussions, A. Biondo, R. H. Brazzle, J. W. 
Epstein, T. Hardt, R. H. Nichols and D. A. Tanner for technical assistance and materials, 
and the Office of Naval Research, DARPA and NSF AST 94-16904 for support. 
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